Social Distancing Correlations Analysis for San Jose
Simone Speizer
4/30/2020
library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(leaflet)
library(lehdr)
options(
tigris_class = "sf",
tigris_use_cache = TRUE
)
Sys.setenv(CENSUS_KEY="10dcd73d7c043e91bac9fb8d3989cbff54b08790")Obtaining demographic variables
Here I obtain various demographic data, including income (percent below 50% and 80% of area median income), vehicle ownership, age, English language ability, and occupants per room.
# obtain the saved census data
setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
acs_vars = readRDS("censusData2018_acs_acs5.rds")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")# load in income data - code adapted from other students
sj_median_income_by_block <-
getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "B19013_001E"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
rename(
Median_Income = B19013_001E
) %>%
filter(!is.na(Median_Income)) %>%
left_join(sj_blockgroups, by = c("blockgroup" = "GEOID")) %>% #this code gives each blockgroup a district designation
filter(
!is.na(DISTRICTS)
) %>%
# this code joins our census data with the social distancing data, processed as shown below
left_join(sj_socialdistancing %>%
filter(weekend == F) %>%
filter(date > shelter_start) %>%
group_by(origin_census_block_group) %>%
summarize(
completely_home_device_count = sum(completely_home_device_count),
device_count = sum(device_count)) %>%
mutate(`% Completely at Home` = (completely_home_device_count/device_count*100) %>% round(1),
`% not completely at home` = (100 - `% Completely at Home`)),
by = c("blockgroup" = "origin_census_block_group")
) %>%
filter(
!is.na(device_count)
) %>%
left_join(sj_pre_sd_at_home_average %>% dplyr::select(origin_census_block_group, `% Completely at Home Pre Shelter`, `% not completely at home pre shelter`), by = c("blockgroup" = "origin_census_block_group"))
sj_ami_by_block <-
getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B19001)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
group_by(blockgroup) %>%
summarize(
Total = B19001_001E,
`Under 75,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E),
#sum(lapply(2:12, function(x) as.name(paste0("B19001_00",x,"E"))))
`Under 100,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E, B19001_013E),
`Under 125,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E, B19001_013E, B19001_014E)
) %>%
mutate(
`% under 75,000` = `Under 75,000` / Total * 100,
`% over 75,000` = (100 - `% under 75,000`),
`% under 100,000` = `Under 100,000` / Total * 100,
`% over 100,000` = (100 - `% under 100,000`),
`% under 125,000` = `Under 125,000` / Total * 100,
`% over 125,000` = (100 - `% under 125,000`),
) %>%
left_join(sj_median_income_by_block %>% dplyr::select(-Median_Income)
) %>%
filter(!is.na(device_count))# loading in language data - code adapted from other students
sj_lang_by_block <-
getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B16004)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(
key = "variable",
value = "estimate",
- blockgroup
) %>%
left_join(acs_vars, by = c("variable" = "name")) %>%
mutate(
tier = substr(label,lapply(label, function(x) max(unlist(gregexpr('!!',x)))+2),nchar(label))
) %>%
filter(tier %in% c('Speak English "not well"',
'Speak English "not at all"',
'Total', 'Speak Spanish',
'Speak Asian and Pacific Island languages')) %>%
group_by(blockgroup, tier) %>%
summarise(
estimate1 = sum(estimate)
) %>%
spread(
key = "tier",
value = "estimate1"
) %>%
mutate(
`% speaking english < well` = (`Speak English "not well"` + `Speak English "not at all"`) / Total * 100,
`% speaking english > well` = (100 - `% speaking english < well`),
`% speaking spanish` = (`Speak Spanish`/ Total) * 100,
`% not speaking spanish` = (100 - `% speaking spanish`),
`% speaking api` = (`Speak Asian and Pacific Island languages` / Total) * 100
) %>%
left_join(sj_median_income_by_block %>% dplyr::select(-Median_Income)) %>%
filter(!is.na(device_count)) %>%
mutate(log_perc = log(`% speaking english < well`))# loading in age data - specifically looking at percentage 65+ and percentage <30
sj_age_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B01001)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(
key = "variable",
value = "estimate",
- blockgroup
) %>%
mutate(
label = acs_vars$label[match(variable,acs_vars$name)]
) %>%
select(-variable) %>%
separate(
label,
into = c(NA,NA,"sex","age"),
sep = "!!"
) %>% filter(!is.na(age)) %>%
mutate(elderly = ifelse(age %in% c("65 and 66 years", "67 to 69 years", "70 to 74 years", "75 to 79 years", "80 to 84 years", "85 years and over"), estimate, NA), `less than 30` = ifelse(age %in% c("Under 5 years", "5 to 9 years", "10 to 14 years", "15 to 17 years", "18 and 19 years", "20 years", "21 years", "22 to 24 years", "25 to 29 years"), estimate, NA)) %>%
group_by(blockgroup) %>%
summarize(elderly = sum(elderly, na.rm = T), `less than 30` = sum(`less than 30`, na.rm = T), total = sum(estimate, na.rm = T)) %>%
mutate(`percent elderly` = elderly*100 / total, `percent less than 30` = `less than 30`*100 / total, `percent nonelderly` = (100 - `percent elderly`)) %>%
left_join(sj_median_income_by_block %>% dplyr::select(-Median_Income)) %>%
filter(!is.na(device_count)) # get data on vehicles available as vehicles allocation
sj_vehicles_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B992512)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
dplyr::select(B992512_001E, blockgroup) %>%
rename(total_vehicles = B992512_001E, blockgroup = blockgroup) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
mutate(`vehicles per capita` = total_vehicles / total) %>%
filter(!is.na(device_count))
# also get data on vehicles available as households without a vehicle
sj_no_vehicles_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B25044)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(key = "variable", value = "estimate", -blockgroup) %>%
mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>%
select(-variable) %>%
separate(label, into = c(NA, NA, NA,"vehicles"), sep = "!!") %>%
filter(!is.na(vehicles)) %>%
group_by(blockgroup, vehicles) %>%
summarize(grouped_vehicles = sum(estimate)) %>%
spread(key = vehicles, value = grouped_vehicles) %>%
mutate(total_nums = `1 vehicle available` + `2 vehicles available` + `3 vehicles available` + `4 vehicles available` + `5 or more vehicles available` + `No vehicle available`, `percent no vehicles` = `No vehicle available`*100 / total_nums, `percent with vehicles` = (100-`percent no vehicles`)) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
filter(!is.na(device_count))# get data on occupants per room
sj_occupants_per_room_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B25014)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(key = "variable", value = "estimate", -blockgroup) %>%
mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>%
select(-variable) %>%
separate(label, into = c(NA, NA, NA,"occupants per room"), sep = "!!") %>%
filter(!is.na(`occupants per room`)) %>%
group_by(blockgroup, `occupants per room`) %>%
summarize(estimate_tot = sum(estimate)) %>%
spread(key = `occupants per room`, value = estimate_tot) %>%
mutate(total_nums = `0.50 or less occupants per room` + `0.51 to 1.00 occupants per room` + `1.01 to 1.50 occupants per room` + `1.51 to 2.00 occupants per room` + `2.01 or more occupants per room`, `percent 1 or more` = (`1.01 to 1.50 occupants per room` + `1.51 to 2.00 occupants per room` + `2.01 or more occupants per room`) * 100/ total_nums, `percent less than 1` = (100-`percent 1 or more`)) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
filter(!is.na(device_count)) Testing correlations
In the plots below, I show the selected variables against percent of devices completely at home since the shelter-in-place order started, as well as against percent of devices pre-shelter-in-place for comparison.
Age:
# age
sj_age_by_block %>%
ggplot(aes(
x = `percent less than 30`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of residents younger than 30",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Young Age Groups"
)
young_model <- lm(sj_age_by_block$`% not completely at home` ~ sj_age_by_block$`percent less than 30`)
summary(young_model)##
## Call:
## lm(formula = sj_age_by_block$`% not completely at home` ~ sj_age_by_block$`percent less than 30`)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.159 -4.833 -0.233 4.473 39.272
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.52788 1.52328 27.919 < 2e-16 ***
## sj_age_by_block$`percent less than 30` 0.20725 0.03861 5.367 1.17e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.352 on 567 degrees of freedom
## Multiple R-squared: 0.04835, Adjusted R-squared: 0.04667
## F-statistic: 28.81 on 1 and 567 DF, p-value: 1.167e-07sj_age_by_block %>% filter(`percent elderly` < 50) %>% # get rid of extreme outliers
ggplot(aes(
x = `percent elderly`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of residents 65 and older",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Elderly Population"
)
elderly_model <- lm(`% not completely at home` ~ `percent elderly`, sj_age_by_block %>% filter(`percent elderly` < 50))
summary(elderly_model)##
## Call:
## lm(formula = `% not completely at home` ~ `percent elderly`,
## data = sj_age_by_block %>% filter(`percent elderly` < 50))
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.119 -5.108 -0.259 4.680 36.632
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.97809 0.79188 66.90 < 2e-16 ***
## `percent elderly` -0.19522 0.05468 -3.57 0.000387 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.456 on 564 degrees of freedom
## Multiple R-squared: 0.0221, Adjusted R-squared: 0.02037
## F-statistic: 12.75 on 1 and 564 DF, p-value: 0.0003869# compare this to pre-shelter-in-place behavior
sj_age_by_block %>%
ggplot(aes(
x = `percent less than 30`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of residents younger than 30",
y = "Percent devices leaving home pre-shelter-in-place",
title = "San Jose: Staying at Home and Young Age Groups Pre Shelter-in-Place"
)
young_model2 <- lm(sj_age_by_block$`% not completely at home pre shelter` ~ sj_age_by_block$`percent less than 30`)
summary(young_model2)##
## Call:
## lm(formula = sj_age_by_block$`% not completely at home pre shelter` ~
## sj_age_by_block$`percent less than 30`)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.0472 -2.4524 -0.1802 2.9456 15.9202
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.50158 0.77360 106.647 < 2e-16 ***
## sj_age_by_block$`percent less than 30` -0.09199 0.01961 -4.691 3.41e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.242 on 567 degrees of freedom
## Multiple R-squared: 0.03736, Adjusted R-squared: 0.03566
## F-statistic: 22.01 on 1 and 567 DF, p-value: 3.41e-06sj_age_by_block %>% filter(`percent elderly` < 50) %>% # get rid of extreme outliers
ggplot(aes(
x = `percent elderly`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of residents 65 and older",
y = "Percent devices leaving home on weekdays pre-shelter-in-place",
title = "San Jose: Staying at Home and Elderly Population Pre Shelter-in-Place"
)
elderly_model2 <- lm(`% not completely at home pre shelter` ~ `percent elderly`, sj_age_by_block %>% filter(`percent elderly` < 50))
summary(elderly_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent elderly`,
## data = sj_age_by_block %>% filter(`percent elderly` < 50))
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.2336 -2.5173 -0.1708 3.0147 12.2215
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.44185 0.39776 194.694 < 2e-16 ***
## `percent elderly` 0.11900 0.02746 4.333 1.74e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.247 on 564 degrees of freedom
## Multiple R-squared: 0.03222, Adjusted R-squared: 0.0305
## F-statistic: 18.77 on 1 and 564 DF, p-value: 1.742e-05Income:
# income - less than $75000
sj_ami_by_block %>%
ggplot(aes(
x = `% over 75,000`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Households Above 50% AMI"
)
income_75_model <- lm(`% not completely at home` ~ `% over 75,000`, sj_ami_by_block)
summary(income_75_model)##
## Call:
## lm(formula = `% not completely at home` ~ `% over 75,000`, data = sj_ami_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.690 -4.649 -0.541 4.168 34.789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.3419 1.1156 57.68 <2e-16 ***
## `% over 75,000` -0.2233 0.0172 -12.99 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.447 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2296, Adjusted R-squared: 0.2282
## F-statistic: 168.7 on 1 and 566 DF, p-value: < 2.2e-16# income - less than $100000
sj_ami_by_block %>%
ggplot(aes(
x = `% over 100,000`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $100,000 (80% AMI) annually",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Households Below 80% AMI"
)
income_100_model <- lm(`% not completely at home` ~ `% over 100,000`, sj_ami_by_block)
summary(income_100_model)##
## Call:
## lm(formula = `% not completely at home` ~ `% over 100,000`, data = sj_ami_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.120 -4.507 -0.736 4.007 33.026
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 61.69920 0.86857 71.03 <2e-16 ***
## `% over 100,000` -0.22077 0.01592 -13.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.33 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2536, Adjusted R-squared: 0.2523
## F-statistic: 192.3 on 1 and 566 DF, p-value: < 2.2e-16# income - less than $125000
sj_ami_by_block %>%
ggplot(aes(
x = `% over 125,000`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $125,000 annually",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Households Below $125,000"
)
income_125_model <- lm(`% not completely at home` ~ `% over 125,000`, sj_ami_by_block)
summary(income_125_model)##
## Call:
## lm(formula = `% not completely at home` ~ `% over 125,000`, data = sj_ami_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.135 -4.477 -0.883 4.279 31.958
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.89204 0.72846 82.22 <2e-16 ***
## `% over 125,000` -0.22989 0.01608 -14.30 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.272 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2653, Adjusted R-squared: 0.264
## F-statistic: 204.4 on 1 and 566 DF, p-value: < 2.2e-16# compare to pre shelter in place
sj_ami_by_block %>%
ggplot(aes(
x = `% over 75,000`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
y = "Percent devices leaving home on weekdays pre-shelter-in-place",
title = "San Jose: Staying at Home and Households Above 50% AMI Pre Shelter-in-Place"
)
income_75_model2 <- lm(`% not completely at home pre shelter` ~ `% over 75,000`, sj_ami_by_block)
summary(income_75_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% over 75,000`,
## data = sj_ami_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.3258 -2.4328 0.0752 2.8199 14.2150
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 74.335377 0.615515 120.769 < 2e-16 ***
## `% over 75,000` 0.074344 0.009487 7.836 2.32e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.109 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.09787, Adjusted R-squared: 0.09627
## F-statistic: 61.4 on 1 and 566 DF, p-value: 2.315e-14# income - less than $100000
sj_ami_by_block %>%
ggplot(aes(
x = `% over 100,000`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $100,000 (80% AMI) annually",
y = "Percent devices leaving home on weekdays pre-shelter-in-place",
title = "San Jose: Staying Home and Households Below 80% AMI Pre Shelter-in-Place"
)
income_100_model2 <- lm(`% not completely at home pre shelter` ~ `% over 100,000`, sj_ami_by_block)
summary(income_100_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% over 100,000`,
## data = sj_ami_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.3698 -2.4317 0.0472 2.7043 14.3737
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.026275 0.481081 155.953 <2e-16 ***
## `% over 100,000` 0.077203 0.008818 8.755 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.06 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1193, Adjusted R-squared: 0.1177
## F-statistic: 76.66 on 1 and 566 DF, p-value: < 2.2e-16# over 125000
sj_ami_by_block %>%
ggplot(aes(
x = `% over 125,000`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $125,000 annually",
y = "Percent devices leaving home on weekdays pre-shelter-in-place",
title = "San Jose: Social Distancing and Households Below $125,000 Pre Shelter-in-Place"
)
income_125_model2 <- lm(`% not completely at home pre shelter` ~ `% over 125,000`, sj_ami_by_block)
summary(income_125_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% over 125,000`,
## data = sj_ami_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.2220 -2.3676 0.1555 2.5383 13.9673
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.432749 0.401316 187.964 <2e-16 ***
## `% over 125,000` 0.085872 0.008858 9.694 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.006 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1424, Adjusted R-squared: 0.1409
## F-statistic: 93.98 on 1 and 566 DF, p-value: < 2.2e-16Language:
# language
sj_lang_by_block %>%
ggplot(aes(
x = `% speaking english > well`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of individuals speaking English well",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and English Language Ability"
)
english_ability_model <- lm(`% not completely at home` ~ `% speaking english > well`, sj_lang_by_block)
summary(english_ability_model)##
## Call:
## lm(formula = `% not completely at home` ~ `% speaking english > well`,
## data = sj_lang_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.526 -4.999 -0.397 3.914 38.388
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.36942 3.36653 20.903 < 2e-16 ***
## `% speaking english > well` -0.22415 0.03775 -5.938 5.02e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.307 on 567 degrees of freedom
## Multiple R-squared: 0.05855, Adjusted R-squared: 0.05689
## F-statistic: 35.26 on 1 and 567 DF, p-value: 5.02e-09sj_lang_by_block %>%
ggplot(aes(
x = `% not speaking spanish`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of individuals not speaking Spanish",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Spanish Language Ability"
)
spanish_speaking_model <- lm(`% not completely at home` ~ `% not speaking spanish`, sj_lang_by_block)
summary(spanish_speaking_model)##
## Call:
## lm(formula = `% not completely at home` ~ `% not speaking spanish`,
## data = sj_lang_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.705 -4.658 -0.741 4.051 37.548
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.80076 1.32061 48.31 <2e-16 ***
## `% not speaking spanish` -0.17128 0.01645 -10.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.845 on 567 degrees of freedom
## Multiple R-squared: 0.1605, Adjusted R-squared: 0.159
## F-statistic: 108.4 on 1 and 567 DF, p-value: < 2.2e-16# compare to pre shelter in place
sj_lang_by_block %>%
ggplot(aes(
x = `% speaking english > well`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of individuals speaking English well",
y = "Percent devices leaving home on weekdays pre-shelter-in-place",
title = "San Jose: Staying at Home and English Language Ability Pre Shelter-in-Place"
)
english_ability_model2 <- lm(`% not completely at home pre shelter` ~ `% speaking english > well`, sj_lang_by_block)
summary(english_ability_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% speaking english > well`,
## data = sj_lang_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.7048 -2.2760 0.0598 2.8560 10.3112
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.51050 1.64215 39.284 <2e-16 ***
## `% speaking english > well` 0.16300 0.01841 8.853 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.052 on 567 degrees of freedom
## Multiple R-squared: 0.1214, Adjusted R-squared: 0.1199
## F-statistic: 78.37 on 1 and 567 DF, p-value: < 2.2e-16sj_lang_by_block %>%
ggplot(aes(
x = `% not speaking spanish`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of individuals not speaking Spanish",
y = "Percent devices leaving home on weekdays pre shelter-in-place",
title = "San Jose: Staying at Home and Spanish Language Ability Pre Shelter-in-Place"
)
spanish_speaking_model2 <- lm(`% not completely at home pre shelter` ~ `% not speaking spanish`, sj_lang_by_block)
summary(spanish_speaking_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% not speaking spanish`,
## data = sj_lang_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.6257 -2.4274 0.0822 2.7537 10.4706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.221669 0.683770 107.08 <2e-16 ***
## `% not speaking spanish` 0.073935 0.008518 8.68 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.062 on 567 degrees of freedom
## Multiple R-squared: 0.1173, Adjusted R-squared: 0.1157
## F-statistic: 75.34 on 1 and 567 DF, p-value: < 2.2e-16Occupants per room:
# occupants per room
sj_occupants_per_room_by_block %>%
ggplot(aes(
x = `percent less than 1`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of households with 1 or fewer occupant per room",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Room Occupancy"
)
occupants_model <- lm(`% not completely at home` ~ `percent less than 1`, sj_occupants_per_room_by_block)
summary(occupants_model)##
## Call:
## lm(formula = `% not completely at home` ~ `percent less than 1`,
## data = sj_occupants_per_room_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.842 -4.849 -0.286 4.137 34.977
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 71.54113 2.97589 24.04 < 2e-16 ***
## `percent less than 1` -0.23399 0.03277 -7.14 2.89e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.126 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.08262, Adjusted R-squared: 0.081
## F-statistic: 50.97 on 1 and 566 DF, p-value: 2.885e-12# compare to pre shelter in place
sj_occupants_per_room_by_block %>%
ggplot(aes(
x = `percent less than 1`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of households with 1 or fewer occupant per room",
y = "Percent devices leaving home on weekdays pre shelter-in-place",
title = "San Jose: Staying at Home and Room Occupancy Pre Shelter-in-Place"
)
occupants_model2 <- lm(`% not completely at home pre shelter` ~ `percent less than 1`, sj_occupants_per_room_by_block)
summary(occupants_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent less than 1`,
## data = sj_occupants_per_room_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.1950 -2.5198 -0.0566 2.8717 14.3938
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65.90546 1.48474 44.389 <2e-16 ***
## `percent less than 1` 0.14478 0.01635 8.854 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.054 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1217, Adjusted R-squared: 0.1201
## F-statistic: 78.4 on 1 and 566 DF, p-value: < 2.2e-16Vehicle ownership:
# vehicles
sj_vehicles_by_block %>%
ggplot(aes(
x = `vehicles per capita`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Number of vehicles per capita",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Vehicles Per Capita"
)
# vehicles - percent with no vehicles
sj_no_vehicles_by_block %>%
ggplot(aes(
x = `percent with vehicles`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with vehicles available",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Vehicle Availability"
)
vehicles_model <- lm(`% not completely at home` ~ `percent with vehicles`, sj_no_vehicles_by_block)
summary(vehicles_model)##
## Call:
## lm(formula = `% not completely at home` ~ `percent with vehicles`,
## data = sj_no_vehicles_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.696 -5.096 -0.235 4.625 38.404
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.75974 5.18753 15.183 < 2e-16 ***
## `percent with vehicles` -0.29763 0.05439 -5.473 6.67e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.268 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.05025, Adjusted R-squared: 0.04858
## F-statistic: 29.95 on 1 and 566 DF, p-value: 6.672e-08# compare to pre shelter in place
sj_no_vehicles_by_block %>%
ggplot(aes(
x = `percent with vehicles`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with vehicles available",
y = "Percent devices leaving home on weekdays pre shelter-in-place",
title = "San Jose: Social Distancing and Vehicle Availability Pre Shelter-in-Place"
)
vehicles_model2 <- lm(`% not completely at home pre shelter` ~ `percent with vehicles`, sj_no_vehicles_by_block)
summary(vehicles_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent with vehicles`,
## data = sj_no_vehicles_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.6695 -2.7766 0.0556 2.9241 10.8094
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65.70738 2.65606 24.739 < 2e-16 ***
## `percent with vehicles` 0.13931 0.02785 5.003 7.55e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.233 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.04235, Adjusted R-squared: 0.04066
## F-statistic: 25.03 on 1 and 566 DF, p-value: 7.55e-07Multiple regression analysis with income, age, language, and occupants per room
# multiple regression
modeltest <- lm(sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` + sj_age_by_block$`percent less than 30` + sj_lang_by_block$`% speaking english > well` + sj_occupants_per_room_by_block$`percent less than 1`)
summary(modeltest)##
## Call:
## lm(formula = sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` +
## sj_age_by_block$`percent less than 30` + sj_lang_by_block$`% speaking english > well` +
## sj_occupants_per_room_by_block$`percent less than 1`)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.884 -4.568 -0.779 4.270 31.930
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 50.772738 4.743391
## sj_ami_by_block$`% over 125,000` -0.250139 0.021391
## sj_age_by_block$`percent less than 30` 0.041453 0.041836
## sj_lang_by_block$`% speaking english > well` 0.090269 0.045832
## sj_occupants_per_room_by_block$`percent less than 1` 0.003916 0.045858
## t value Pr(>|t|)
## (Intercept) 10.704 <2e-16 ***
## sj_ami_by_block$`% over 125,000` -11.694 <2e-16 ***
## sj_age_by_block$`percent less than 30` 0.991 0.3222
## sj_lang_by_block$`% speaking english > well` 1.970 0.0494 *
## sj_occupants_per_room_by_block$`percent less than 1` 0.085 0.9320
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.253 on 563 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.273, Adjusted R-squared: 0.2679
## F-statistic: 52.86 on 4 and 563 DF, p-value: < 2.2e-16Additions based on literature search
I also consider education and internet access, based on previous research. Education:
sj_education_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B15003)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(key = "variable", value = "estimate", -blockgroup) %>%
mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>%
select(-variable) %>%
separate(label, into = c(NA, NA, "education level"), sep = "!!") %>%
mutate(`education level` = replace_na(`education level`, "total_educ")) %>% # if the education level field is NA, this corresponded to the total number in that blockgroup
spread(key = `education level`, value = estimate) %>%
mutate(`percent associates or higher` = (`Associate's degree` + `Bachelor's degree` + `Doctorate degree` + `Master's degree`)*100/total_educ, `percent less than associates` = 100-`percent associates or higher`) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
filter(!is.na(device_count))
# plotting
sj_education_by_block %>%
ggplot(aes(
x = `percent associates or higher`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of people with an degree at Associate's level or higher",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Education"
)
educ_model <- lm(`% not completely at home` ~ `percent associates or higher`, sj_education_by_block)
summary(educ_model)##
## Call:
## lm(formula = `% not completely at home` ~ `percent associates or higher`,
## data = sj_education_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.749 -4.612 -0.937 3.826 42.245
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.23755 0.83655 72.01 <2e-16 ***
## `percent associates or higher` -0.20682 0.01642 -12.60 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.568 on 567 degrees of freedom
## Multiple R-squared: 0.2187, Adjusted R-squared: 0.2174
## F-statistic: 158.7 on 1 and 567 DF, p-value: < 2.2e-16# compare to pre shelter in place
sj_education_by_block %>%
ggplot(aes(
x = `percent associates or higher`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of people with an degree at Associate's level or higher",
y = "Percent devices leaving home on weekdays pre-shelter-in-place",
title = "San Jose: Social Distancing and Education Pre Shelter-in-Place"
)
educ_model2 <- lm(`% not completely at home pre shelter` ~ `percent associates or higher`, sj_education_by_block)
summary(educ_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent associates or higher`,
## data = sj_education_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.8071 -2.2967 0.1566 2.6958 11.5796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.470385 0.450731 167.440 < 2e-16 ***
## `percent associates or higher` 0.074211 0.008844 8.391 3.86e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.078 on 567 degrees of freedom
## Multiple R-squared: 0.1105, Adjusted R-squared: 0.1089
## F-statistic: 70.4 on 1 and 567 DF, p-value: 3.863e-16Internet:
sj_internet_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B28002)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(key = "variable", value = "estimate", -blockgroup) %>%
mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>%
select(-variable) %>%
separate(label, into = c(NA, NA, "subscription", "type", "additional"), sep = "!!") %>%
filter(is.na(subscription) | (type == "Broadband such as cable, fiber optic or DSL") & is.na(additional)) %>%
mutate(type = replace_na(type, "total_num")) %>%
dplyr::select(blockgroup, type, estimate) %>%
spread(key = type, value = estimate) %>%
mutate(`percent high speed` = `Broadband such as cable, fiber optic or DSL`*100/total_num, `percent no high speed` = 100-`percent high speed`) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
filter(!is.na(device_count))
# plotting
sj_internet_by_block %>%
ggplot(aes(
x = `percent high speed`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of households with broadband such as cable, fiber optic or DSL",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and High Speed Internet"
)
internet_model <- lm(`% not completely at home` ~ `percent high speed`, sj_internet_by_block)
summary(internet_model)##
## Call:
## lm(formula = `% not completely at home` ~ `percent high speed`,
## data = sj_internet_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.988 -4.550 -0.406 3.971 38.924
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.06603 2.24315 32.57 <2e-16 ***
## `percent high speed` -0.28052 0.02751 -10.20 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.798 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1552, Adjusted R-squared: 0.1537
## F-statistic: 104 on 1 and 566 DF, p-value: < 2.2e-16# compare to pre-shelter-in-place behavior
sj_internet_by_block %>%
ggplot(aes(
x = `percent high speed`,
y = `% not completely at home pre shelter`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of households with broadband such as cable, fiber optic or DSL",
y = "Percent devices leaving home on weekdays pre-shelter-in-place",
title = "San Jose: Social Distancing and High Speed Internet Pre Shelter-in-Place"
)
internet_model2 <- lm(`% not completely at home pre shelter` ~ `percent high speed`, sj_internet_by_block)
summary(internet_model2)##
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent high speed`,
## data = sj_internet_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.4202 -2.6369 -0.0996 2.7840 14.1874
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 71.52221 1.20356 59.425 < 2e-16 ***
## `percent high speed` 0.09226 0.01476 6.251 8.02e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.184 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.06458, Adjusted R-squared: 0.06293
## F-statistic: 39.08 on 1 and 566 DF, p-value: 8.023e-10Do another multiple regression analysis, this time with education and income
educ_income_model <- lm(sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` + sj_education_by_block$`percent associates or higher`)
summary(educ_income_model)##
## Call:
## lm(formula = sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` +
## sj_education_by_block$`percent associates or higher`)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.770 -4.458 -0.815 4.137 33.078
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 61.46824 0.81750
## sj_ami_by_block$`% over 125,000` -0.16632 0.02233
## sj_education_by_block$`percent associates or higher` -0.08880 0.02194
## t value Pr(>|t|)
## (Intercept) 75.191 < 2e-16 ***
## sj_ami_by_block$`% over 125,000` -7.449 3.53e-13 ***
## sj_education_by_block$`percent associates or higher` -4.047 5.91e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.175 on 565 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.286, Adjusted R-squared: 0.2835
## F-statistic: 113.2 on 2 and 565 DF, p-value: < 2.2e-16Try with internet and income
educ_income_model <- lm(sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` + sj_internet_by_block$`percent high speed`)
summary(educ_income_model)##
## Call:
## lm(formula = sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` +
## sj_internet_by_block$`percent high speed`)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.317 -4.498 -0.731 4.156 31.808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.55452 2.26972 28.442 <2e-16
## sj_ami_by_block$`% over 125,000` -0.20017 0.02109 -9.491 <2e-16
## sj_internet_by_block$`percent high speed` -0.07295 0.03365 -2.168 0.0306
##
## (Intercept) ***
## sj_ami_by_block$`% over 125,000` ***
## sj_internet_by_block$`percent high speed` *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.248 on 565 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2714, Adjusted R-squared: 0.2688
## F-statistic: 105.2 on 2 and 565 DF, p-value: < 2.2e-16Testing animating the plot
# collect the demographic variables
sj_dem_distancing <- sj_internet_by_block %>%
dplyr::select(`percent high speed`, `% not completely at home`, blockgroup) %>%
left_join(sj_education_by_block %>% dplyr::select(blockgroup, `percent associates or higher`)) %>%
left_join(sj_ami_by_block %>% dplyr::select(blockgroup, `% over 125,000`)) %>%
left_join(sj_age_by_block %>% dplyr::select(blockgroup, `percent less than 30`)) %>%
left_join(sj_lang_by_block %>% dplyr::select(blockgroup, `% not speaking spanish`)) %>%
left_join(sj_no_vehicles_by_block %>% dplyr::select(blockgroup, `percent with vehicles`)) %>%
left_join(sj_occupants_per_room_by_block %>% dplyr::select(blockgroup, `percent less than 1`))
# another collection for pre shelter in place behavior
sj_dem_distancing_pre_shelter <- sj_dem_distancing %>%
dplyr::select(-`% not completely at home`) %>%
left_join(sj_internet_by_block %>% dplyr::select(`% not completely at home pre shelter`, blockgroup))
# relabel column for leaving home
colnames(sj_dem_distancing_pre_shelter)[9] <- "% not completely at home"
sj_dem_distancing[is.na(sj_dem_distancing)] <- 0
sj_dem_distancing_pre_shelter[is.na(sj_dem_distancing_pre_shelter)] <- 0
# add column indicating before or after shelter in place, then bind the two sets of data
sj_dem_distancing_pre_shelter <- sj_dem_distancing_pre_shelter %>%
mutate(
income_trendline =
fitted(lm(sj_dem_distancing_pre_shelter$`% not completely at home` ~ sj_dem_distancing_pre_shelter$`% over 125,000`))
) %>%
cbind(`Before or After Shelter-in-Place` = "before")
sj_dem_distancing <-
sj_dem_distancing %>%
mutate(
income_trendline =
fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`% over 125,000`))
) %>%
cbind(`Before or After Shelter-in-Place` = "after") %>%
rbind(sj_dem_distancing_pre_shelter)
# try animating
fig <-
plot_ly (sj_dem_distancing) %>%
add_trace(
x = ~`% over 125,000`,
y = ~`% not completely at home`,
frame = ~`Before or After Shelter-in-Place`,
type = 'scatter',
mode = 'markers'
) %>%
add_trace(
x = ~`% over 125,000`,
y = ~income_trendline,
type = 'scatter',
mode = 'lines',
line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
frame = ~`Before or After Shelter-in-Place`
) %>%
animation_button(visible = F)
fig# # save as rds
# saveRDS(sj_dem_distancing, "/Users/simonespeizer/pCloud Drive/Shared/SFBI/Restricted Data Library/Safegraph/covid19analysis/sj_sd_dem_data.rds")
# fig <- plot_ly(sj_dem_distancing) %>%
# add_trace(
# x = ~`% over 125,000`,
# y = ~`% not completely at home`,
# frame = ~`Before or After Shelter-in-Place`,
# type = "scatter",
# mode = "markers",
# name = "Under $125,000",
# marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
# visible = T
# ) %>%
# add_trace(
# x = ~`% over 125,000`,
# y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`% over 125,000`)),
# name = 'trendline',
# mode = 'lines',
# line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
# frame = ~`Before or After Shelter-in-Place`,
# visible = F
# ) %>%
# add_trace(
# x = ~`% not speaking spanish`,
# y = ~`% not completely at home`,
# frame = ~`Before or After Shelter-in-Place`,
# name = "speak spanish",
# marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
# visible = F
# ) %>%
# add_trace(
# x = ~`% not speaking spanish`,
# y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`% not speaking spanish`)),
# name = 'trendline',
# mode = 'lines',
# line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
# frame = ~`Before or After Shelter-in-Place`,
# visible = F
# ) %>%
# add_trace(
# x = ~`percent associates or higher`,
# y = ~`% not completely at home`,
# frame = ~`Before or After Shelter-in-Place`,
# name = "percent higher degree",
# marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
# visible = F
# ) %>%
# add_trace(
# x = ~`percent associates or higher`,
# y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`percent associates or higher`)),
# name = 'trendline',
# mode = 'lines',
# line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
# frame = ~`Before or After Shelter-in-Place`,
# visible = F
# ) %>%
# add_trace(
# x = ~`percent high speed`,
# y = ~`% not completely at home`,
# frame = ~`Before or After Shelter-in-Place`,
# name = "percent high speed internet access",
# marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
# visible = F
# ) %>%
# add_trace(
# x = ~`percent high speed`,
# y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`percent high speed`)),
# name = 'trendline',
# mode = 'lines',
# line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
# frame = ~sj_dem_distancing$`Before or After Shelter-in-Place`,
# visible = F
# ) %>%
# add_trace(
# x = ~`percent less than 30`,
# y = ~`% not completely at home`,
# frame = ~`Before or After Shelter-in-Place`,
# name = "percent less than 30",
# marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
# visible = F
# ) %>%
# add_trace(
# x = ~`percent less than 30`,
# y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`percent less than 30`)),
# name = 'trendline',
# mode = 'lines',
# line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
# frame = ~`Before or After Shelter-in-Place`,
# visible = F
# ) %>%
# layout(
# updatemenus = list(
# list(
# active = 0,
# type = 'buttons',
# buttons = list(
# list(
# label = "Households Under $125,000",
# method = 'update',
# args = list(list(visible = c(T, T, F, F, F, F, F, F, F, F)),
# list(title = "Under $125,000",
# xaxis = list(title = "% Households Under $125,000 in Income")))),
# list(
# label = "Speaking Spanish",
# method = 'update',
# args = list(list(visible = c(F, F, T, T, F, F, F, F, F, F)),
# list(title = "Not Speaking Spanish",
# xaxis = list(title = "% Residents Not Speaking Spanish")))),
# list(
# label = "Education Level",
# method = 'update',
# args= list(list(visible = c(F, F, F, F, T, T, F, F, F, F)),
# list(xaxis = list(title = "% Residents With Associate's Degree or Higher")))),
# list(
# label = "High Speed Internet",
# method = 'update',
# args= list(list(visible = c(F, F, F, F, F, F, T, T, F, F)),
# list(xaxis = list(title = "% Households With High Speed Internet Access")))),
# list(
# label = "Young Population",
# method = 'update',
# args= list(list(visible = c(F, F, F, F, F, F, T, T, F, F)),
# list(xaxis = list(title = "% Residents Under Age 30"))))
# )
# )
# ),
# yaxis = list(title = "% Residents Leaving Home",
# font = list(size = 15)),
# showlegend = FALSE
# )
# figExperimentation
Experimentation with other variables and other ways of analyzing the social distancing data. First I look at a few other possible variables. I start with units in the structure.
# try getting other variables
# get data on units in structure
sj_units_in_structure_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B25024)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(key = "variable", value = "estimate", -blockgroup) %>%
mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>%
select(-variable) %>%
separate(label, into = c(NA, NA, "units"), sep = "!!") %>%
filter(!is.na(units)) %>%
spread(key = units, value = estimate) %>%
mutate(total_nums = `1, attached` + `1, detached` + `10 to 19` + `2` + `20 to 49`+ `3 or 4` + `5 to 9`+ `50 or more`+ `Boat, RV, van, etc.`+ `Mobile home`, `percent 20 or more` = (`20 to 49`+`50 or more`)* 100/ total_nums, `percent 1 only` = (`1, attached` + `1, detached`)*100/total_nums, `percent > 1` = 100 - `percent 1 only`) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
filter(!is.na(device_count))
# plot
sj_units_in_structure_by_block %>%
ggplot(aes(
x = `percent 20 or more`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of structures with more than 20 units",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and 20 or More Units Per Structure"
)
summary(lm(`% not completely at home` ~ `percent 20 or more`, sj_units_in_structure_by_block))##
## Call:
## lm(formula = `% not completely at home` ~ `percent 20 or more`,
## data = sj_units_in_structure_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.765 -5.190 -0.201 4.843 37.335
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.06531 0.40940 122.290 <2e-16 ***
## `percent 20 or more` 0.03712 0.02052 1.809 0.071 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.459 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.005748, Adjusted R-squared: 0.003992
## F-statistic: 3.272 on 1 and 566 DF, p-value: 0.07099sj_units_in_structure_by_block %>%
ggplot(aes(
x = `percent 1 only`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of structures with only one unit",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Only 1 Unit Per Structure"
)
summary(lm(`% not completely at home` ~ `percent 1 only`, sj_units_in_structure_by_block))##
## Call:
## lm(formula = `% not completely at home` ~ `percent 1 only`, data = sj_units_in_structure_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.608 -5.078 -0.224 4.405 38.261
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.24549 0.88920 62.130 < 2e-16 ***
## `percent 1 only` -0.06648 0.01132 -5.872 7.33e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.237 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.05743, Adjusted R-squared: 0.05576
## F-statistic: 34.48 on 1 and 566 DF, p-value: 7.328e-09Household type and size:
# load data on household type and size
sj_house_size_type_by_block <- getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = "group(B11016)"
) %>%
mutate(
blockgroup =
paste0(state,county,tract,block_group)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(key = "variable", value = "estimate", -blockgroup) %>%
mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>%
select(-variable) %>%
separate(label, into = c(NA, NA, "type", "size"), sep = "!!") %>%
filter(!is.na(type))
# household type
sj_house_type_by_block <- sj_house_size_type_by_block %>%
filter(is.na(size)) %>%
dplyr::select(-size) %>%
spread(key = type, value = estimate) %>%
mutate(`total households` = `Family households` + `Nonfamily households`, `percent nonfamily` = `Nonfamily households` / `total households`) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
filter(!is.na(device_count))
sj_house_type_by_block %>%
ggplot(aes(
x = `percent nonfamily`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent nonfamily households",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Household Type"
)
summary(lm(`% not completely at home` ~ `percent nonfamily`, sj_house_type_by_block))##
## Call:
## lm(formula = `% not completely at home` ~ `percent nonfamily`,
## data = sj_house_type_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.828 -5.089 0.052 4.597 38.208
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.7692 0.6403 74.609 < 2e-16 ***
## `percent nonfamily` 11.0270 2.2223 4.962 9.24e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.305 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.04169, Adjusted R-squared: 0.03999
## F-statistic: 24.62 on 1 and 566 DF, p-value: 9.245e-07# household size
sj_house_size_by_block <- sj_house_size_type_by_block %>%
filter(!is.na(size)) %>%
dplyr::select(-type) %>%
group_by(blockgroup, size) %>%
summarize(`total of this size` = sum(estimate)) %>%
spread(key = size, value = `total of this size`) %>%
mutate(total_nums = `1-person household` + `2-person household` + `3-person household` + `4-person household` + `5-person household`+ `6-person household` + `7-or-more person household`, `percent 5 or more` = (`5-person household`+`6-person household` + `7-or-more person household`)* 100/ total_nums, `percent 1 or 2 only` = (`1-person household` + `2-person household`)*100/total_nums) %>%
left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
filter(!is.na(device_count))
sj_house_size_by_block %>%
ggplot(aes(
x = `percent 5 or more`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of households with 5 or more people",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Households With 5 or More"
)
summary(lm(`% not completely at home` ~ `percent 5 or more`, sj_house_size_by_block))##
## Call:
## lm(formula = `% not completely at home` ~ `percent 5 or more`,
## data = sj_house_size_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.419 -4.958 -0.605 4.419 37.588
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.67686 0.56617 85.976 < 2e-16 ***
## `percent 5 or more` 0.10054 0.02541 3.957 8.55e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.369 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.02692, Adjusted R-squared: 0.0252
## F-statistic: 15.66 on 1 and 566 DF, p-value: 8.545e-05sj_house_size_by_block %>%
ggplot(aes(
x = `percent 1 or 2 only`,
y = `% not completely at home`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of households with 1 or 2 people",
y = "Percent devices leaving home on weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Small Household Size"
)
summary(lm(`% not completely at home` ~ `percent 1 or 2 only`, sj_house_size_by_block))##
## Call:
## lm(formula = `% not completely at home` ~ `percent 1 or 2 only`,
## data = sj_house_size_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.408 -5.383 -0.223 4.811 37.673
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.45494 0.98251 50.335 <2e-16 ***
## `percent 1 or 2 only` 0.02185 0.02043 1.069 0.285
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.475 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.002016, Adjusted R-squared: 0.0002529
## F-statistic: 1.143 on 1 and 566 DF, p-value: 0.2854Next I consider different ways of looking at the social distancing data. First I try distance traveled.
# try other ways of looking at the social distancing data
# first look at total distance traveled
sj_sd_distance <- sj_socialdistancing %>%
filter(date > shelter_start) %>%
group_by(origin_census_block_group) %>%
summarize(total_dist_traveled = sum(distance_traveled_from_home), device_count = sum(device_count)) %>%
mutate(total_dist_per_device = total_dist_traveled / device_count)
sj_distance_testing <- left_join(sj_ami_by_block, sj_sd_distance, by = c("blockgroup" = "origin_census_block_group")) %>% left_join(sj_age_by_block %>% select(blockgroup, `percent less than 30`))
sj_distance_testing %>% filter(total_dist_per_device < 500) %>%
ggplot(aes(
x = `% over 75,000`,
y = total_dist_per_device
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
y = "Average distance traveled per device during weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Income, Distance Metric"
)
This is very skewed by outliers, and probably not a useful metric.
Now I consider including devices that traveled <1km as staying at (or near) home.
sj_sd_range <- sj_socialdistancing %>%
filter(weekend == F) %>%
filter(date > shelter_start) %>%
mutate(travel_buckets_split = lapply(bucketed_distance_traveled, function(x) strsplit(x, "<1000")[[1]][2]), less_than_1km = lapply(travel_buckets_split, function(x) strsplit(x, ":")[[1]][2]), less_than_1km = lapply(less_than_1km, function(x) strsplit(x, ",")[[1]][1])) %>%
mutate(less_than_1km = lapply(less_than_1km, function(x) str_remove(x, "[}]"))) %>% # clean a bit more
mutate(less_than_1km = as.numeric(less_than_1km), less_than_1km = replace_na(less_than_1km, 0)) %>%
mutate(home_or_1km = completely_home_device_count + less_than_1km) %>%
group_by(origin_census_block_group) %>%
summarize(home_or_1km = sum(home_or_1km), device_count = sum(device_count)) %>%
mutate(`% Within 1km of Home` = (home_or_1km/device_count*100) %>% round(1), `% farther than 1km` = (100-`% Within 1km of Home`))
# join this with other data
sj_1km_testing <- left_join(sj_ami_by_block, sj_sd_range, by = c("blockgroup" = "origin_census_block_group")) %>%
left_join(sj_occupants_per_room_by_block %>% dplyr::select(`percent less than 1`, blockgroup)) %>%
left_join(sj_age_by_block %>% dplyr::select(`percent less than 30`, blockgroup)) %>%
left_join(sj_lang_by_block %>% dplyr::select(`% speaking english > well`, blockgroup))
# plot with income
sj_1km_testing %>%
ggplot(aes(
x = `% over 75,000`,
y = `% farther than 1km`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
y = "% of devices going farther than 1km of home, weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Income, 1km Range"
)
summary(lm(`% farther than 1km` ~ `% over 75,000`, sj_1km_testing))##
## Call:
## lm(formula = `% farther than 1km` ~ `% over 75,000`, data = sj_1km_testing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.159 -4.993 -0.401 4.266 41.458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.2525 1.1743 47.90 <2e-16 ***
## `% over 75,000` -0.2343 0.0181 -12.95 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.839 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2285, Adjusted R-squared: 0.2271
## F-statistic: 167.6 on 1 and 566 DF, p-value: < 2.2e-16# plot with age
sj_1km_testing %>%
ggplot(aes(
x = `percent less than 30`,
y = `% farther than 1km`
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of people younger than 30",
y = "Percent of devices farther than 1km of home during weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Age, 1km Range"
)
summary(lm(`% farther than 1km` ~ `percent less than 30`, sj_1km_testing))##
## Call:
## lm(formula = `% farther than 1km` ~ `percent less than 30`, data = sj_1km_testing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.413 -5.309 -0.304 4.766 42.271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.21800 1.57767 19.787 < 2e-16 ***
## `percent less than 30` 0.27326 0.03999 6.833 2.15e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.651 on 567 degrees of freedom
## Multiple R-squared: 0.07608, Adjusted R-squared: 0.07445
## F-statistic: 46.69 on 1 and 567 DF, p-value: 2.154e-11# run multiple regression model
modeltest2 <- lm(sj_1km_testing$`% farther than 1km` ~ sj_1km_testing$`% over 75,000` + sj_1km_testing$`percent less than 30` + sj_1km_testing$`% speaking english > well` + sj_1km_testing$`percent less than 1`)
summary(modeltest2)##
## Call:
## lm(formula = sj_1km_testing$`% farther than 1km` ~ sj_1km_testing$`% over 75,000` +
## sj_1km_testing$`percent less than 30` + sj_1km_testing$`% speaking english > well` +
## sj_1km_testing$`percent less than 1`)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.400 -4.765 -0.793 4.599 40.972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.12037 4.91129 10.205 < 2e-16
## sj_1km_testing$`% over 75,000` -0.20967 0.02267 -9.249 < 2e-16
## sj_1km_testing$`percent less than 30` 0.14057 0.04466 3.148 0.00173
## sj_1km_testing$`% speaking english > well` -0.02227 0.04822 -0.462 0.64432
## sj_1km_testing$`percent less than 1` 0.01303 0.04897 0.266 0.79024
##
## (Intercept) ***
## sj_1km_testing$`% over 75,000` ***
## sj_1km_testing$`percent less than 30` **
## sj_1km_testing$`% speaking english > well`
## sj_1km_testing$`percent less than 1`
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.77 on 563 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.246, Adjusted R-squared: 0.2406
## F-statistic: 45.92 on 4 and 563 DF, p-value: < 2.2e-16It looks like the fit of these selected variables is slightly better for the social distancing data based on not traveling farther than 1km.
Now I also consider “non-work” behavior.
sj_nonworking_by_block <- sj_socialdistancing %>%
filter(weekend == F) %>%
filter(date > shelter_start) %>%
mutate(nonworking = device_count - completely_home_device_count - part_time_work_behavior_devices - full_time_work_behavior_devices) %>%
group_by(origin_census_block_group) %>%
summarize(nonworking_count = sum(nonworking), total_device = sum(device_count)) %>%
mutate(nonworking_percent = nonworking_count*100 / total_device, percent_only_work_home = 100-nonworking_percent) %>%
left_join(sj_1km_testing %>% dplyr::select(`% over 75,000`, `percent less than 30`, `% speaking english > well`, `percent less than 1`, blockgroup), by = c("origin_census_block_group" = "blockgroup"))
# plot against age and income
sj_nonworking_by_block %>%
ggplot(aes(
x = `% over 75,000`,
y = nonworking_percent
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
y = "Percent of devices leaving home for non-work purposes during weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Income, Nonworking Behavior"
)
summary(lm(nonworking_percent ~ `% over 75,000`, sj_nonworking_by_block))##
## Call:
## lm(formula = nonworking_percent ~ `% over 75,000`, data = sj_nonworking_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.1765 -3.1343 -0.1581 3.1028 17.2184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.24081 0.77032 53.54 <2e-16 ***
## `% over 75,000` -0.12775 0.01187 -10.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.142 on 566 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1698, Adjusted R-squared: 0.1683
## F-statistic: 115.8 on 1 and 566 DF, p-value: < 2.2e-16sj_nonworking_by_block %>%
ggplot(aes(
x = `percent less than 30`,
y = nonworking_percent
)) + geom_point() +
geom_smooth(method=lm) +
labs(
x = "Percent of people younger than 30",
y = "Percent of devices leaving home for non-work purposes during weekdays since shelter-in-place",
title = "San Jose: Social Distancing and Age, Nonworking Behavior"
)
summary(lm(nonworking_percent ~ `percent less than 30`, sj_nonworking_by_block))##
## Call:
## lm(formula = nonworking_percent ~ `percent less than 30`, data = sj_nonworking_by_block)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.6311 -3.5775 -0.2204 3.4844 17.7223
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.16288 1.02739 30.332 <2e-16 ***
## `percent less than 30` 0.05479 0.02604 2.104 0.0358 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.633 on 567 degrees of freedom
## Multiple R-squared: 0.007747, Adjusted R-squared: 0.005997
## F-statistic: 4.427 on 1 and 567 DF, p-value: 0.03582# multiple regression model
modeltest3 <- lm(sj_nonworking_by_block$nonworking_percent ~ sj_nonworking_by_block$`% over 75,000` + sj_nonworking_by_block$`percent less than 30` + sj_nonworking_by_block$`% speaking english > well` + sj_nonworking_by_block$`percent less than 1`)
summary(modeltest3)##
## Call:
## lm(formula = sj_nonworking_by_block$nonworking_percent ~ sj_nonworking_by_block$`% over 75,000` +
## sj_nonworking_by_block$`percent less than 30` + sj_nonworking_by_block$`% speaking english > well` +
## sj_nonworking_by_block$`percent less than 1`)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.9975 -3.2297 0.0295 3.1285 16.9204
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 52.91517 3.21912 16.438
## sj_nonworking_by_block$`% over 75,000` -0.10453 0.01486 -7.035
## sj_nonworking_by_block$`percent less than 30` -0.07166 0.02927 -2.448
## sj_nonworking_by_block$`% speaking english > well` -0.04189 0.03161 -1.326
## sj_nonworking_by_block$`percent less than 1` -0.07375 0.03210 -2.298
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## sj_nonworking_by_block$`% over 75,000` 5.82e-12 ***
## sj_nonworking_by_block$`percent less than 30` 0.0147 *
## sj_nonworking_by_block$`% speaking english > well` 0.1855
## sj_nonworking_by_block$`percent less than 1` 0.0219 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.093 on 563 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.19, Adjusted R-squared: 0.1842
## F-statistic: 33.01 on 4 and 563 DF, p-value: < 2.2e-16These variables do worse for the percent nonworking metric, which makes sense.